Assessment of a Hybrid Method for Hypersonic Flows

A preliminary hybrid particle-continuum computational framework for simulating hypersonic interacting flows is proposed. The framework consists of the direct simulation Monte Carlo-Information Preservation (DSMC-IP) method coupled with a Navier-Stokes solver. Since the DSMC-IP method provides the macroscopic information in each time step, determination of the continuum fluxes across the interface between the particle and continuum domains becomes straightforward. A hypersonic flow over a two-dimensional wedge is considered as an example and compared with pure particle calculations. The results show that this preliminary hybrid framework is promising but several issues are yet to be resolved. INTRODUCTION In recent years, biconic configurations and hollow cylinder/flares were extensively studied with different numerical approaches [1, 2, 3, 4] for a set of hypersonic free-stream condition. In comparisons with experimental data [5], it was evident that the computational fluid dynamics (CFD) techniques based on the Navier-Stokes (NS) equations have better performance than the direct simulation Monte Carlo (DSMC) method [6] in terms of the capture of flow structures and the prediction of the separation region size. However, it was shown that the CFD methods constantly over predict the heat transfer on the fore body ahead of where separation occurs. Furthermore, it was demonstrated [7] that adjacent to the body wall, especially near the cone tip, there is a strong thermal nonequilibrium effect, that is inconsistent with the basic assumptions of the NS equations. In addition, Roy et al. [8] showed that DSMC and CFD do agree for conditions that are more rarefied than the experiments. Therefore, the objective of the present work is to develop a combined DSMC-CFD computational method for the physically accurate and numerically efficient analysis of this class of hypersonic flows. Two primary issues must be taken into account in developing a combination of DSMC and CFD methods. We first need to know when to switch between the methods. Since it is well known that the NS equations are not valid under rarefied conditions, it is general to use a continuum breakdown parameter as the criterion. For the hypersonic flows mentioned above, this issue has been investigated in a previous study [9] and where it was concluded that a proposed parameter Knmax max KnD KnT KnV (1) with a threshold value of 0.05 can best predict the regions where the Navier-Stokes equations fail. The Knudsen number in Eq. 1 is expressed as KnQ λ Q ∇Q , with Q being any flow property in general. Since continuum breakdown is often related to the transport phenomena of viscosity and heat transfer, we only consider the flow properties of density(D), temperature(T) and magnitude of velocity(V). Another issue regarding the development of a hybrid method is the approach for information exchange at the interface between the DSMC and CFD domains. At the interface, macroscopic flow properties must be provided to the CFD method to evaluate the net fluxes and to the DSMC method to initialize the particles entering from the continuum domain into the rarefaction domain. Several approaches have been considered, such as the Marshak condition [10], the kinetic flux-vector splitting (KFVS) scheme [11, 12] and the adaptive mesh and algorithm refinement (AMAR) embedding a particle method [13]. Unfortunately, the development of a robust, multi-dimensional scheme that is capable of handling nonequilibrium, hypersonic compressed flows has not yet been accomplished. The primary difficulty of the second issue results from the fact that the DSMC approach always involves very strong statistical scatter unless the number of samples is large enough. Recently, the information preservation method (IP) was proposed to reduce the statistical scatter in low-speed, constant density flow systems [14]. Since then, the method has been generalized to allow density and temperature variations and shown to be effective in solving microfluidic and micro-electro-mechanical system (MEMS) problems [15, 16, 17]. One advantage of the IP method is that the macroscopic values of the flow field are known at any time, since the information in cells is updated for each time step. This eases the complexity of coupling the particle method with the continuum solver. A hybrid approach that combines the IP method and a NS solver has shown great progress in solving micro-scale gas flows [18]. This motivates us in the present study to explore the IP method for hypersonic, interacting flows. The layout of the paper is as follows. A brief description of the continuum approach will be presented first in the next section, followed by an introduction of the DSMC-IP method. A detailed explanation on how to combine these two methods together is provided in the section of Domain Coupling. A hypersonic flow over a two-dimensional wedge is considered to assess the new hybrid technique. In the last section, conclusions and suggestions for future work are provided. NUMERICAL SCHEMES Continuum Approach The Navier-Stokes equations in the continuum domain are solved numerically using an explicit Gauss-Seidel line relaxation method and second-order accurate, modified Steger-Warming flux vector splitting [19]. The viscosity μ is modeled with the power law and the thermal conductivity κ is determined from the Prandtl number μ μref T Tref ω Pr cpμ κ where cp is the specific heat at constant pressure. A slip-boundary model proposed by Gökçen [20] is implemented. DSMC-IP Approach The information preservation method was first developed by Fan and Shen [14] to overcome the statistical scatter problem in DSMC simulations, especially for systems in which the flow speed is much smaller than the molecular speed. In addition to the ordinary thermal velocity that is utilized to compute the particle trajectory, each simulation particle in the DSMC-IP method also possesses macroscopic preserved information such as velocity vector and temperature. The DSMC-IP method has achieved great success for solving micro-scale gas flows [see 15, 16, 17]. In the most recent work by Sun and Boyd [21], an additional temperature term is introduced in order to solve the contradiction between the real flux and the DSMC-IP representation of the translational energy flux across a cell interface. From gas kinetic theory, the average translational energy of a molecule at equilibrium temperature T is 3kT 2 (k is the Boltzmann constant), whereas the average translational energy carried by a molecule across an interface is 2kT . The extra energy must be correctly modeled or the energy of the whole system will not be conserved. Therefore, the simulation particles in the IP method have an additional temperature Ta. In each time step of the DSMC-IP method, simulation particles are first moved and collided in the usual way as in the standard DSMC method. The preserved velocity in the ri direction and temperature of simulation particles are updated by solving ∂Vi ∂ t 1 ρc ∂ pc ∂ ri (2) ∂ ∂ t V 2 i 2 ξ R T 2 1 ρc ∂ ∂ ri Vi c pc (3) where p is the pressure, ρ is the mass density, ξ is the number of internal degrees of freedom of molecules, and the subscript c denotes the macroscopic information for the computational cells. After the preserved information of simulation particles is updated, the preserved information for cells is updated by taking the arithmetic average over the information of all Np particles in the cell. Vi c 1 Np Np ∑ j 1 Vi j (4) Tc 1 Np Np ∑ j 1 Tj Ta j (5) The density is updated by solving the continuity equation ∂ρc ∂ t ∂ ∂ ri ρc Vi c (6) The ideal gas law, p ρRT , is assumed. A detailed description and implementation of the DSMC-IP method can be found in Ref. [21]. Note that an adequate numerical scheme must be employed to solve Eq. 6 due to the presence of shock waves in supersonic flows. Since the continuity equation also appears in the NS equations, we use the same technique described in the last subsection for solving it. The current DSMC-IP code is based on a parallel optimized DSMC code named MONACO [22]. A sub-cell scheme is implemented for selection of collision pairs where the number of sub-cells is scaled by the local mean free path. DOMAIN COUPLING To implement the coupling between the particle method and the NS solver, buffer and reservoir DSMC-IP cells are introduced in the continuum domain adjacent to the domain interface, as depicted in Fig. 1. A similar concept of reservoir cells was first proposed in Ref. [23]. The buffer DSMC-IP cells work as an extension of the particle domain. Simulation particles that end their movement phase within the pure particle domain or in the buffer cells are retained. Those that leave these two regions are removed. For each time step, all simulation particles in the reservoir cells are first deleted and then re-generated based on the cell-centered values. The number of new particles is evaluated from the cell density value and the particle velocities and temperature are initialized to the Chapman-Enskog distribution [24] based on the corresponding cell values. The newly generated particles are randomly distributed within the reservoir cells. In this study, one layer of buffer cells and five layers of reservoir cells are employed. In the continuum domain, the NS solver determines the interface continuum fluxes by using the NS variables and DSMC-IP cell macroscopic information. Since the macroscopic information in the DSMC-IP cells is known in each time step, the DSMC-IP cells adjacent to the domain interface are treated as the ghost cells that provide the boundary conditions for flux computations. NUMERICAL EXAMPLE This section describes a Mach 4 numerical experiment to assess the hybrid technique. A geometry of a 2D wedge with a 25 half angle is chosen. The symmetric line of the wedge is aligned with the free-stream. Therefore, only the upper half of the wedge is considered. Since the wake region behind the wedge is not of interest in this investigation, we assume that the wedge is infinitely long but only the first 5 cm from the leading edge is considered. A structured grid, 300 cells along the body by 200 cells normal to the body, used in all our computations is shown in Fig. 2. The fluid is pure nitrogen and the free-stream conditions are: U∞ 1111 1 m/s, T∞ 185 6 K and ρ∞ 6 545 10 4 kg/m3. The mean-free-path in the free-stream is about 10 4 m. An isothermal wall at a temperature of 293.3 K is assumed. These specific flow conditions are the same as those in CUBRC Run 28 [25] except for a lower flow speed. In the CFD calculation, μref 1 656 10 5 N s/m2, Tref 273 K and ω 0 74. The Prandtl number is considered as a constant of 0.72. A pure DSMC steady state solution is obtained with the use of more than 1.8 million simulation particles at the end of the computation. The reference time step in the pure DSMC calculation is 5 nsec. 300,000 time steps of computation are performed and the last 50,000 time steps are sampled to obtain the results. In the hybrid simulation, a steady state solution from the CFD method is first obtained. Using the steady state solution and Knmax, one can determine the locations of the interfaces between the continuum and the particle domains. In the particle region, cell values are set to the CFD steady state results and simulation particles are initialized to the Chapman-Enskog distribution. Continuum Domain Particle Domain Interface Buffer Cells Reservoir Cells Interface Pure NS Cells Pure DSMC-IP Cells Extended DSMC-IP Cells NS Step DSMC-IP Step Ghost-Cells from DSMC-IP FIGURE 1. Interface cell types. X (cm) Y (c m ) 0 1 2 3 4 5 0 1 2 3 4 5 300x200 Cells (every 5th point shown in each direction) FIGURE 2. Grid employed for 25 wedge.